library('SDMTools')

mc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/'
bc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/BC6_monthly/'
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/microCLIM_images/'

mcfiles = list.files(mc.dir,full.names=T, pattern='.asc.gz')
mcnames = gsub('_','',gsub('.asc.gz','',list.files(mc.dir, pattern='.asc.gz')))

bcfiles = list.files(bc.dir,full.names=T, pattern='.asc.gz')[1:4]
bcnames = gsub('.asc.gz','',list.files(bc.dir,pattern='.asc.gz')[1:4])

for (i in 1:length(mcfiles))

	{
	
	### Read in two ASCII's and find their residual
	
	t.mc = read.asc.gz(mcfiles[i])
	
	t.bc = read.asc.gz(bcfiles[i])
	
	t.diff = t.mc - t.bc
	
	### Establish break points for color scheme NB needs to be checked when looping, these breaks ARE NOT equally spaced.... i.e. won't work in filled.contour command
	
	t.breaks = seq(min(t.diff,na.rm=T),0,(max(t.diff,na.rm=T)-min(t.diff,na.rm=T))/100)
	t.breaks2 = seq(0,max(t.diff,na.rm=T),(max(t.diff,na.rm=T)-min(t.diff,na.rm=T))/100)
	t.breaks2 = seq(0,max(t.diff,na.rm=T),(max(t.diff,na.rm=T)-min(t.diff,na.rm=T))/100)[2:length(t.breaks2)]
	zero.calc = abs(max(t.breaks)-min(t.breaks2))/3
	breaks.fin = c(t.breaks,max(t.breaks)+zero.calc,min(t.breaks2)-zero.calc,t.breaks2)
	
	#### Establish color palettes
	
	cols.below=colorRampPalette(c('purple','blue','green','white'))
	cols.above = colorRampPalette(c('white','yellow','orange','red'))
	
	### Get a list of colors that corresponds to the breaks
	
	cb.codes = cols.below(length(which(breaks.fin<0)))
	ca.codes = cols.above(length(which(breaks.fin>0)))[2:length(which(breaks.fin>0))]
	col.codes = c(cb.codes,ca.codes)

	
	png(paste(out.dir,'MC_minus_BC Surface_0',substr(bcnames[i],4,4),'.png',sep=''),units='in',width=8,height=11,res=500)
	
	image(t.diff, breaks = breaks.fin, col=col.codes,oldstyle=T) # NB breaks must be one element longer than cols
	
	pnts=cbind(x=c(147,147.25,147.25,147),y=c(-16,-16,-17,-17))
	
	legend.gradient(pnts,cols=col.codes,limits = c(round(min(t.diff,na.rm=T),2),round(max(t.diff,na.rm=T),2)),title=paste('Surface 0',substr(bcnames[i],4,4),sep=''))
	
	dev.off()
	
	cat(paste('\n',bcnames[i],' - Image Completed','\n',sep=''))
	
	}
	
# Done

col.levs.below = seq(round(min(zz$z,na.rm=T)-1,0),0,.25)+.125
col.levs.above = seq(0,round(max(zz$z,na.rm=T)+1,0),.25)-.125
levs = c(col.levs.below,col.levs.above[3:length(col.levs.above)])

### Create a color pallete

cols.below = colorRampPalette(c('blue','green','white'))
cols.above = colorRampPalette(c('white','orange','red'))

col.codes.below = cols.below(length(col.levs.below)-1)
col.codes.above = cols.above(length(col.levs.above[3:length(col.levs.above)]))



### Now plot the object zz using the image command

png(paste(out.dir,'persptest.png',sep=''))
image(zz)
filled.contour(zz,levels=levs,col = c(col.codes.below,col.codes.above),xlab="dem",ylab="fpc")
dev.off()



